Critical fronts in initiation of excitation waves 
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We consider the problem of initiation of propagating wave in a one- dimensional excitable fiber. 
In the FitzHugh-Nagumo theory, the key role is played by "critical nucleus" and "critical pulse" 
solutions whose (center-)stable manifold is the threshold surface separating initial conditions leading 
to propagation and those leading to decay. We present evidence that in cardiac excitation models, 
this role is played by "critical front" solutions. 



PACS numbers: 87.19.Hh, 87.19.La, 02.90.+P 



I. INTRODUCTION 

An excitable medium is a thermodynamically non- 
equilibrium system that has a stable spatially uniform 
"resting state", but responds to an above-threshold lo- 
calized stimulus by a propagating non-decaying "exci- 
tation wave" . Excitation waves play key roles in living 
organisms and are observed in chemical and physical sys- 
tems, e.g. nerves, heart muscle, catalytic redox reactions, 
large aspect lasers and star formation in galaxies [l[ . Un- 
derstanding conditions of successful initiation of excita- 
tion waves is particularly important for heart where such 
waves trigger coordinated contraction of the muscle and 
where a failure of initiation can cause or contribute to 
serious or fatal medical conditions, or render inefficient 
the work of pacemakers or defibrillators [|J . 

The theoretical understanding of excitability stems 
from FitzHugh's simplified model of a nerve mem- 
brane 3 ■ One of his key concepts is "quasi-threshold" , 
which gets precise in the limit of large time scale sepa- 
ration between the processes of excitation and recovery. 
Then the fast subsystem has unstable "threshold" equi- 
libria. Initial conditions below such an equilibrium lead 
to decay, and those above it lead to excitation. 

In a spatially extended FitzHugh-Nagumo (FHN) sys- 
tem 0, [i[ , the ability of a stimulus to initiate a wave 
depends also on its spatial extent, the aspect summa- 
rized by Rushton's [|| concept of "liminal length" . More 
generic is the concept of "critical curve" in the stimulus 
strength-spatial extent plane (see Fig. [T]) . A stimulus ini- 
tiates a wave if its parameters are above this curve, or 
leads to decay if below. 

Mathematically, the problem is about classification of 
initial conditions that will or will not lead to a travel- 
ing wave solution. The key question is the nature of the 
boundary between the two classes. A detailed analysis of 
this boundary has been done for the FitzHugh-Nagumo 
system and its variations. This has led to the concept of a 
critical nucleus, discussed below in more detail. Roughly, 
this is a spatially extended analog of a threshold equilib- 
rium in the point system: critical nucleus is also a sta- 
tionary but unstable solution, and its small perturbations 
lead to either initiation of an excitation wave, for pertur- 
bations of one direction, or to decay, for perturbations of 



the opposite direction. 

We stress that although the role of FHN as a universal 
prototype of excitable systems has been disputed, to our 
knowledge, there are still no alternatives to the critical 
nucleus concept, as far as initiation problem is concerned. 

In this paper, we present evidence that cardiac excita- 
tion provides an example of alternative type of system, in 
which there is no place for the critical nucleus. Two inde- 
pendent observations led to this study. Firstly, numerical 
simulations of the cardiac excitation models reveal signif- 
icant qualitative differences in the way initiation occurs 
in such models, compared to the FHN-style systems [f|. 
Secondly, asymptotic analysis of detailed cardiac exci- 
tation models reveals that in the fast subsystem there, 
there is no analoqoi the unstable threshold equilibrium 
of FHN systems [7j , and the threshold there has a com- 
pletely different mathematical nature. Further, elemen- 
tary arguments show that in cardiac equations there are 
no nontrivial stationary solutions that could play the role 
similar to the critical nucleus in FHN system. 

Thus we have a theoretical vacuum here. Obviously, 
one cannot even begin to think about investigating ini- 
tiation criteria without understanding the nature of the 
critical solutions. This paper aims to fill this vacuum, 
and clarify the nature of the critical solutions. We an- 
alyze a simplified model of cardiac excitation, and use 
the knowledge of its exact solutions to demonstrate that 
for this model the concept of critical nucleus should be 
replaced with a new concept of critical front. We also 
confirm numerically the relevance of this new concept on 
an example of a detailed ionic cardiac excitation model. 



II. FITZHUGH-NAGUMO SYSTEM 

First we recapitulate some known theoretical concepts 
related to initiation of waves (see e.g. H,[| and references 
therein). We consider FitzHugh-Nagumo system in the 
form 



u xx + f(u) 
e(au — v) 



/(u) = u(u-0)(l-u) 



(1) 



where e > 0, a > 0, 9 6 (0, 1/2) (some works consider 
piecewise linear functions / of similar shape) on a half- 
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fiber, (x, t) G [0, oo) x [0, oo) with a no-flux boundary 

Ux(0,t)=0 > (2) 

and a rectangular initial perturbation of width a; st i m and 
amplitude M s tim, 



u{x, 0) 



i9(a: B 



:), t;(a:,0) = (3) 



where O(-) is the Heaviside step function. 

Fig. [2a,b) shows two typical results of the initiation 
process: a successful initiation, leading to generation of a 
propagating pulse, and an unsuccessful, leading to decay 
of excitation in the whole half-fiber into the resting state. 

If e — 0, problem ()1|2I3|) reduces to an initiation 
problem for Z eldovich- Frank- K amene t sky ( Z FK ) equa- 
tion also known as Nagumo equation [ill ]. 



u t = u xx + f(u), 
u x (0,t) = 0, 

u(x,0) = u stim Q(x s ti m - x). 



(4) 
(5) 
(6) 



Fig. [ljc,d) illustrates the initiation and its failure in this 
reduced problem. Instead of a propagating pulse, suc- 
cessful initiation produces a propagating front. In the full 
model with small e, this front is followed, in time scale 
0£ _1 , by a wave-back to form a full excitation pulse. 

A key role in understanding initiation belongs to 
an unstable nontrivial bounded time-independent so- 
lution of (|4I5[) . sometimes called critical nucleus, 
by analogy with phase transitions theory. Such 
solution is unique; for a cubical nonlinearity / 
as in (fT]), this solution has the form u CT (x) = 



36\/2 (l + 6>)\/2 + cosh(2;Ve)\/2-56l + 6> 2 . Its lin- 
earization spectrum has exactly one unstable eigenvalue, 
while all other eigenvalues are stable. So the stable man- 
ifold of this stationary solution has codimension one, and 
divides the phase space of (|4l5p to two open sets. One of 
these sets corresponds to initial conditions leading to suc- 
cessful initiation, and the other to decay. In particular, if 
the initial condition satisfies u(x, 0) < u cr (x), x G [0, oo) 
then u(x,t) decays as t — > oo, and if u{x,Q) > u CI (x), 
x S [0, oo) then u(x, t) approaches a stable propagating 
front solution. Moreover, if a continuous one-paramctric 
family of initial conditions contains some that initiate a 
wave and some that lead to decay, then there is always at 
least one that does neither, but gives a solution that ap- 
proaches the critical nucleus. This critical nucleus is the 
same for all such families, e.g. docs not depend on the 
shape of the initial distribution u(x, 0), as long as its am- 
plitude is at the threshold corresponding to that shape. 
Initial conditions very close to the threshold generate so- 
lutions which approach the critical nucleus and then de- 
part from it, either toward propagation or toward decay. 
This transient stationary state can be seen in Fig. [TJc,d) 
where the initial conditions are selected very close to the 
threshold. 

For small e > 0, system (jTJ does not have nontrivial 
stationary solutions. However, for x € (— oo, oo) equation 



(fT]) has an unstable propagating pulse solution u cr (x— ct), 
v C r(z — ct) such that u C i{x) — > it cr (x), v cr (x) — > and 
c = 0c 1 / 2 as £ \ 0. Due to translational symmetry, this 
solution (in the comoving frame of reference) has a zero 
eigenvalue corresponding to the eigenfunction (u' CI , v' CI ). 
This solution also has a single unstable eigenvalue. So 
its center-stable manifold has codimension one and is the 
threshold hypcrsurface dividing the phase space into the 
decay domain and the initiation domain. So here we have 
a critical pulse solution, which we define as an unsta- 
ble traveling wave that is asymptotic to the resting state 
for both limits x — ct — > ±oo. For small e, the critical 
pulse is essentially a slowly traveling variant of the crit- 
ical nucleus. Any solution with the initial condition at 
the threshold hypcrsurface will asymptotically approach 
this critical pulse (suitably shifted), and any solutions 
starting close to the threshold will approach this critical 
pulse as a transient [2(j. This is illustrated in Fig. [51 

With this understanding, the excitation condition in 
terms of (x s tim> u a ti m ) reduces to computing the intersec- 
tion of the two-parametric manifold described by ^ with 
the codimension 1 stable (center-stable) manifold of the 
critical nucleus (critical pulse). This gives the curve on 
the (Xstimj itgtim) plane separating initial conditions lead- 
ing to excitation propagation from those leading to decay. 
This can be done numerically or, with appropriate sim- 
plifications, analytically. An example of dealing with this 
problem in the ZFK equation, using Galerkin style ap- 
proximations can be found in [l2T |. Here we concentrate 
on the principal question, about the nature of the thresh- 
old hypersurface in the functional space, for cardiac ex- 
citation equations. It appears that in cardiac equations, 
this nature is different from the FitzHugh-Nagumo the- 
ory just considered. 



III. SIMPLIFIED CARDIAC EXCITATION 
MODEL 

Now we consider the simplified model of /Na-driven ex- 
citation fronts in typical cardiac excitation models pro- 
posed in [l3l |: 

E t =E xx + e(E-l)h, h t = (e{-E)-h)/r (7) 
with boundary condition 

E x (0,t) = (8) 

and initial condition 



E(x,0) = -a + E atim G(x a 



x), h(x,0) = l, (9) 



where the variable E represents transmembrane potential 
of the cardiac tissue, h is the probability of the Na-gates 
being open, r is a dimensionless parameter and a > 
represents the pre- frontal voltage which we consider fixed 
in this paper. System can be obtained by simplifying 
right-hand sides of the fast subsystem in an appropriate 
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FIG. 1: Initiation of excitation in FitzHugh-Nagumo system. (a,b) Full system (I I 121.' H "FHN" for parameter values: a — 0.37, 
9 — 0.13, e — 0.02. Stimulation parameters: £ s tim — 2.10 for both, below-threshold u s tim = 0.43, leading to decay, for (a) 
and above-threshold u s tim = 0.44, leading to initiation of excitation propagation, for (b). (c,d) Fast subsystem (|4l5l6p "ZFK": 
same parameters as in (a,b) except e — 0. Stimulation: a; s tim = 2.10 for both, below-threshold Ustim = 0.3304831 for (c) and 
above-threshold u s ti m = 0.3304833 for (d). Bold black lines: initial conditions, (e) The corresponding critical curves, separating 
initiation initial conditions from decay initial conditions. Simulation done on an interval x € [0, L], L = 120 with Neuman 
boundaries, central space differencing with h x — 0.15, and explicit Euler timestepping with h t = 0.01. 



* + lOOu 
500 I i 




(a) (b) (c) (d) 

FIG. 2: The critical pulse is a universal transient for any near-threshold initial condition. The solutions to {TJ for slightly 
below-threshold (a,c) and slightly above-threshold (b,d) amplitudes, for smaller stimulus width x B ti m = 2.10 in (a,b) and 
larger x s t\m = 10.05 in (c,d). Parameter values: e = 0.02, a — 0.37, ht = 0.01, h x = 0.15, L = 120. Stimulus amplitudes: 
,i stim = 0.431929399574766 for (a), 0.431929399574768 for (b), 0.191802079312694 for (c) and 0.191802079312696 for (d). In 
all cases we see a slow, low-amplitude unstable propagating pulse which subsequently either decays or evolves into a fast, 
high-amplitude stable propagating pulse. 



asymptotic limit of a typical cardiac excitation model 
[14| . In that sense, system ([7]) plays the same role for 
a typical cardiac excitation model, as ZFK equation (H|) 
plays for a classical activator-inhibitor excitable system 
like ©. 

System J7]) does not have nontrivial bounded station- 
ary solutions: if E t = h t = then any bounded solution 
has the form E = a, h = 0(— a), a = const. So, there 
are no critical nuclei in this system. Nevertheless, sys- 
tem ([7]) gives propagating front solutions for initial con- 
ditions above a threshold and decay for those below it. 
Hence there is a question, what happens when the initial 
condition is exactly at the threshold. 

System ((7J) has a family of propagating front solutions 



E(z) 
h(z) 




(10) 



where 



x - ct, u = 1 + rc 2 (a + 1), A 



In I 



and parameters c, a and r related by 

tc 2 In ((1 + a){c 2 + t- 1 )) + In (l + a" 1 ) = 0. (11) 

For a fixed a, there is a r*(a) such that for r > r*, 
equation (TIT]) has two solutions for c, c = c±(a, r), c+ > 
c_ [l3[. There is numerical and analytical evidence that 
solutions with c = c + arc stable and those with c = c_ 
are unstable with one positive eigenvalue [l3l . [l5| . 

Hence by analogy with the FHN system, we propose 
the following 

Conjecture 1 The center-stable manifold of the unsta- 
ble front solution mU\) with c = c_(a, r) is the threshold 
hypersurface, separating the initial conditions leadin g to 
initiation from the initial conditions leading to decay, \2j\1 

That is, instead of a critical nucleus or a critical pulse 
solution, the role of the threshold solution is played by 
a "critical front" , which we define as a traveling wave 
solution with different asymptotics at x — ct — > +oo and 
x — ct — > — oo: the pre- frontal state and the post-frontal 
state. 

An "experimentally testable" consequence of this con- 
jecture is that for any initial condition exactly at the 
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threshold, the solution will approach the unstable front 
as t — > +00. For any initial condition near the threshold, 
the solution will come close to the unstable front and stay 
in its vicinity for a long time: if the positive eigenvalue 
is A and the initial condition is <5-close to the threshold, 
the transient front should be observed for the time of 
the order of A - 1 1 In <5| . This transient front solution will 
not depend on the initial condition, as long as the initial 
condition is at the threshold. 



t + loos 
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FIG. 3: Evolution of two different near-threshold initial con- 
ditions toward the critical front solution in system (0. Ini- 
tial stimuli: x sti m = 0.3, E Btim = 12.716330706144868 (upper 
row) andi stim = 1.5, E stim = 2.619968799545055 (lower row). 
Other parameters: r = 8.2, a = 1, h x = 0.075, h t = 0.0025, 
L = 50. 

We have tested these predictions by numerical simula- 
tion of (|7I8I9|) . The results arc shown in Figs. [3] and [4] 

Fig.[3]illustrates two solutions starting from initial con- 
ditions with different x s ti m values. In both cases, -E s tim 
values have been chosen close to the respective thresholds 
with high precision. In both cases, the solutions evolve 
in the long run toward the same propagating front. 

Fig. d] presents an analysis of a pair of solutions, one 
with slightly above-threshold and the other with slightly 
below-thrcshold initial conditions. To separate the evolu- 
tion of the front shape from its movement, we employed 
the idea of symmetry group decomposition with explicit 
representation of the orbit manifold (see e.g. [l6|). Prac- 
tically, we define the front point Xf = Xf(t) via 

E(x f (t),t) = E* 

for some constant E* which is guaranteed to be repre- 
sented exactly once in the front at every instant of time 
(we have chosen E* = 0). Then E(x — Xf(t),t) gives 
the voltage profile "in the standard position" , and x / (t) 
describes the movement of this profile. 

The predictions based on the Conjecture are that the 
voltage profile should, after an initial transient depending 
on the initial condition, approach the profile of the slow 
unstable front solution given by (|10[) with c = c_ (t, a) 
and stay close to it for some time, before either develop- 
ing into the fast stable front (JTUJ) with c = c + (r,a) or 




100 200 300 -50 -25 25 50 200 400 600 

X x - Xf t 

FIG. 4: Transient critical fronts are close to the unstable 
front solution of ([7j). Initial conditions: = 1.5, with 

E s tim = 2.619968799545055 in the upper row and E st i m = 
2.619968799545054 in the lower row, other parameters the 
same as in Fig. [3] Left column: evolution of the E pro- 
files in the laboratory frame of reference. Middle column: 
same evolution, in the frame of reference comoving with the 
front. Right column: speed of the front. Blue/green dashed 
lines in the middle and right columns correspond to the exact 
fast /slow front solutions of ([7j). 



decaying. Likewise, the speed of the front should, after 
an initial transient, be close to the speed of the slow un- 
stable front c_ (r, a) , before either switching the speed of 
the fast stable front c+(r, a) or dropping to zero. This 
is precisely what is seen on Fig. HI where we have taken 
advantage of knowing the exact solutions E(x — c±t) and 
c± for both the fast and the slow fronts. 

Initial conditions with different x st i m and -E s tim close 
to the corresponding threshold produce the same picture 
with the exception of the initial transient. We have also 
checked that the length of the time period during which 
the solution stays close to the unstable front is, roughly, 
a linear function of the number of correct decimal figures 
in -Estim, as it should be according to the Conjecture. 



IV. DETAILED CARDIAC EXCITATION 
MODEL 

The simplified model ([7]) is quantitatively very far from 
any realistic ionic model of cardiac excitation, and has 
many peculiar qualitative features, stemming from the 
non-standard asymptotic embedding leading to it. Hence 
the newly described phenomenon of critical front could 
be an artifact of the simplifications. To eliminate this 
possibility, we have tested the relevance of the critical 
front concept to a full ionic model of cardiac excitation. 
We have chosen the model of human atrial tissue due to 
Courtemanche et al. (CRN) [13], which is less stiff than 
a typical ventricular or Purkinje fiber model, is well for- 
mulated in a mathematical sense and is popular among 
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served critical front may well be the front of a critical 
pulse solution in the full model. However Fig. [5] demon- 
strates that the critical front is a practical and well work- 
ing concept even for the full model, unlike the critical 
pulse which may be theoretically existing but practically 
unobservable: notice the number of significant decimal 
digits in initial conditions required to produce only the 
critical front observed for 80 ms, and recall that the num- 
ber of decimals is roughly proportional to the duration 
of the observation of an unstable solution. 



FIG. 5: Critical fronts in CRN model [17[. Shown are volt- 
age profiles in every 10 ms. Parameter values: h t — 0.01ms, 
h x = 0.2, L = 40, the length unit chosen so that voltage 
diffusion coefficient equals 1. Stimulus witdh s: s ti m = 2, stim- 
ulus amplitudes: V stim = 29.31542299307152 mV (left panel) 
and V s tim = 29.31542299307153 mV (right panel). The crit- 
ical fronts are formed within first 10 ms and then are seen 
for subsequent 80 ms on both panels before exploding into an 
excitation wave of much bigger amplitude and speed on the 
right panel, and decaying on the left panel. 



cardiac modellers. This model operates with 21 dynamic 
variables including the transmembrane voltage V. We 
have used the default parameter values as described in 
fll\ and supplemented the equation for V in the ODE 
system with a diffusion term Dd 2 V/dx 2 . As the spatial 
scale is not important for the question at hand, we as- 
sumed D = 1. The initial conditions for V were taken in 
the form 



V(x,0)=V r + V gtim Q(x B 



x) 



where V r = — 81.18 mV is the standard resting potential, 
and for all other 20 variables at their resting values as 
described in [Tt} ■ Fig. [5] illustrates a pair of solutions 
with initial conditions slightly above and slightly below 
the threshold. The critical front solution is clearly seen 
there: it has the upper voltage of about —46 mV and 
during 80 ms of its existence propagates with a speed 
approximately 0.06 space units per millisecond. Then for 
the above-critical case it develops into an excitation front 
with maximal voltage about +3mV and speed 0.8 space 
units per millisecond, and decays for the below-critical 
case. 

Mathematically, the post-front voltage of about 
—46 mV observed in Fig. [5] is not a true equilibrium of 
the full CRN model so the critical front can only be an 
asymptotic concept in an appropriate asymptotic embed- 
ding, say as ones described in 14 1 or and the ob- 



V. DISCUSSION 

We have presented numerical evidence that the center- 
stable manifold of the unstable slow front solution of the 
fast subsystem of a cardiac excitation model serves as the 
threshold hypcrsurface separating initial conditions lead- 
ing to successful initiation and those leading to decay. 
This means e.g. that a critical curve for a two-parametric 
family of initial conditions, be that family ((9|) with pa- 
rameters (:c s tim, -Estim) or an Y other, can be found as an 
intersection of this codimension-1 critical hyper-surface 
with the two-dimensional manifold of those initial con- 
ditions. Finding this hypersurface can be done numer- 
ically or analytically using suitable approximation, e.g. 
as it was done in [T3, EH for the ZFK equation; this is a 
subject for further investigation. 

Another problem for future study is to verify that the 
findings remain qualitatively true for formal asymptotic 
embeddings of various cardiac excitation models, and in 
particular, which mathematical features of these embed- 
dings are essential for the existence of the critical fronts. 
Of principal importance is the conclusion that for car- 
diac equations, instead of the "critical nucleus" or its 
slowly moving variant "critical pulse" known from the 
FitzHugh-Nagumo theory, we now have a "critical front 
solution" . This, in particular, means physically that the 
make-or-break conditions of cardiac excitation wave are 
restricted to a vicinity of its front. 

Finally, we believe that this study sets a useful ex- 
ample for initiation problems in other types of excitable 
systems, alternative to the existing critical nucleus the- 
ory, since not all, if any, real-world excitable systems are 
well described by an asymptotic structure as in (TT]). 
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